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ABSTRACT 


The  governing  equations  for  laminar,  transitional 
and/or  turbulent,  incompressible,  three-dimensional  boundary 
layers  are  solved  numerically.  The  equations  are  developed 
in  terms  of  an  orthogonal  surface  coordinate  system  and  include 
surface  curvature  effects.  The  coordinates  are  generated  nu- 
merically, while  the  inviscid  flow  is  obtained  from  a general 
pof'^ti^  flow  procedure.  The  only  restrictions  on  body  geometry 
possess  a blunt  nose  and  a plane  of  svmmetry.  The 

bound  ary-^l 


ayer  equations,  after  being  transformed  into  similarity 
tjipe  variables,  are  solved  using  the  implicit  finite-difference 
Krause  scheme.  Various  test  cases  are  presented  to  establish 
the  accuracy  of  the  resulting  computer  program. 
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SECTION  I 


I 


INTRODUCTION 

A computer  program  has  been  developed  to  predict 
laminar  or  turbulent  three-dimensional  boundary  layers  over 
blunt  bodies  inmersed  in  incompressible  fluids.  The  program 
includes  the  effects  of  surface  curvature.  Such  a code  is 
useful  for  predicting  the  boundary- layer  growth  over  sub- 
merged bodies  at  angle  of  attack,  or  bodies  of  non-circular 
cross-section  at  zero  angle  of  attack. 

In  this  report,  the  three-dimensional  boundary- layer 
equations  are  developed  in  a surface -or iented , orthogonal, 
curvilinear  coordinate  system  with  all  surface  curvature  ef- 
fects included.  Further,  the  limiting  forms  of  these  equa- 
tions in  S5anmetry  planes  and  at  the  stagnation  point  are 
developed.  It  is  required  that  any  body  shape  considered  by 
the  code  possess  a plane  of  symmetry. 

The  resulting  equations  are  integrated  using  a marching, 
implicit  finite-difference  scheme  on  an  IBM  370  system-model 
158  digital  computer. 

1 . Background 

The  three-dimensional  boundary- layer  equations  for 
laminar,  compressible  flow  have  been  developed  by  Moore  (Ref.  1) . 
Procedures  for  numerically  solving  these  equations  for  the  in- 
compressible case  have  been  developed  by  Blottner  and  Ellis 
(Ref.  2)  and  Chang  and  Patel  (Ref.  3).  The  Blottner  and  Ellis 
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procedure  for  laminar  three-dimensional  boundary  layers  uses 
an  orthogonal,  surface -oriented,  curvilinear  coordinate  system 
which  is  numerically  generated.  In  this  procedure,  the  co- 
ordinate system  is  developed  independently  of  the  boundary- 
layer  or  inviscid  flow,  but  has  as  its  origin  the  inviscid 
stagnation  point.  The  procedure  does  require,  however,  an 
analytic  solution  for  the  inviscid  flow  and  is  therefore 
restricted  to  a rather  narrow  class  of  bodies.  The  method  of 
Chang  and  Patel  for  laminar  or  turbulent  three-dimensional 
boundary  layers  also  uses  orthogonal  surface  coordinates,  but 
these  coordinates  are  calculated  analytically  based  upon  the 
geometry  of  mathematically  well  defined  shapes.  As  such,  this 
procedure  is  also  limited  in  the  class  of  bodies  which  can  be 
calculated. 

The  procedure  developed  in  the  present  report  takes 
advantage  of  the  versatility  of  the  Blottner-Ellis  (Ref.  2) 
coordinate  system  so  that  a much  wider  class  of  bodies  may  be 
studied.  Essentially,  the  numerically  generated  coordinate 
package  of  the  Blottner-Ellis  procedure  is  separated  from  the 
boundary-layer  procedure.  This  allows  for  the  development  of 
a coordinate  system  for  a given  geometry  to  be  carried  out  as 
a separate  task  from  the  boundary- layer  calculation.  Thus,  once 
the  coordinate  system  has  been  developed  for  a given  geometry 
at  a given  angle  of  attack,  that  body  can  be  studied  under  a 
variety  of  flow  conditions  without  having  to  regenerate  the 
body  coordinate  system.  Further,  since  the  procedure  is 
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nijmerical,  it  can  be  used  with  a discrete  description  of  body 
geometry  as  well  as  analytic,  thus  further  enhancing  the  pro- 
cedures versatility. 

Since  the  body  geometry  that  can  be  handled  by  the 
code  is  restricted  only  to  blunt  shapes  with  a plane  of  S3mi- 
metry,  a rather  general  technique  for  predicting  the  inviscid 
flow  about  the  body  is  required.  The  Smith-Hess  (Ref.  4)  pro- 
cedure is  used  for  this  purpose.  An  intermediate  code  is  then 
used  to  interpolate  the  output  of  the  Smith-Hess  code  onto  the 
surface-orthogonal  coordinate  system.  All  of  these  data  are 
then  Fourier  fitted  and  the  Fourier  coefficients  are  stored  on 
disk.  This  disk  is  then  read  by  the  bound ary -layer  program 
which  predicts  the  boundary-layer  development  on  the  body. 

The  technique  described  in  this  report  accounts  for  all 

I surface  curvature  effects,  and,  to  the  authors'  knowledge,  this 

f 

[ is  the  first  time  these  effects  have  been  accounted  for  in 

i 

three-dimensional  boundary  layers. 
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SECTION  II 


GOVERNING  EQUATIONS 

1.  Coordinate  System  and  Boundary-Layer  Equations 

In  this  section,  the  full  incompressible  three- 
dimensional  boundary -layer  equations  are  developed.  The  equa- 
tions developed  are  written  in  an  orthogonal  curvilinear  body- 
oriented  coordinate  system  whose  origin  is  located  at  the  stag- 
nation point  for  blunt  nosed  bodies.  The  details  of  the  co- 
ordinate system  are  depicted  in  Fig.  1 . As  can  be  seen  in  this 
figure,  two  coordinate  systems  are  shown;  the  first  a body - 
oriented  polar  coordinate  system  and  the  second  a surface  - 
oriented  curvilinear  system.  In  this  system,  C*  Is  measured 
away  from  the  stagnation  point  along  meridional  cuts,  u)*  is 
measured  normal  to  ^and  around  the  body  and  n*  is  measured 
normal  to  the  body.  At  the  stagnation  point  0,  and  a)*=  0 
along  the  windward  symmetry  cut.  In  this  system  the  body  radius 
r*  is  generally  a function  of  w*  For  this  coordinate  system  a 
differential  distance  ds*  is 

ds*^  = h*^  (C*.u)*,n*)d?*^  + h*^(5*,a)*,n*)da)*^  + dn*^  (1) 

where  h*  and  h*  are  the  metric  scale  coefficients. 

C to 

The  governing  equations  of  fluid  mechanics  written  in 

this  coordinate  system  are  (Ref. 5 ): 

Incompressibility : 
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We  DOW  approach  the  above  equations  with  the  assumption  i 

{ 

that  in  a vanishingly  small  region  near  the  surface,  a region  | 

of  order  e thick,  all  flow  properties  except  the  normal  ve-  j 

locity  component  are  of  order  one  ~ v itself  being  of  order  t. 

I 

Thus,  we  introduce  the  following  stretched  coordinates: 

C = I,  0)  * u,  n * n/e  (9) 

where  e « Re”^^^.  According  to  our  above  described  ordering  ' 

scheme  we  then  write  i 

u(r,w,n;e)  » u(|,w,n',e)  ! 

w(5,t»r,n;e)  = w(^,a>,n;e)  I 

vU,w,n:e)  » e~^v(C,  w.nje)  ^ 

p(C,w,n;e)  = p(C,w,n*E)  (10) 


and  note  that  h^  “ h^  and  h^^  * h^,. 

Introduction  of  Eq . (9)  and  (10)  into  the  governing 
equations  and  retention  of  all  terms  to  order  e gives 
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Equations  (11)  through  (15)  represent  the  full  second -order 
three-dimensional  incompressible  boundary-layer  equations 
written  in  the  surface -oriented  curvilinear  coordinate  system 
described  above. 

2.  Boundary  and  Matching  Conditions 

In  order  to  complete  the  above  described  set  of  equa- 
tions, the  boundary  conditions  at  the  body  surface  and  matching 
conditions  at  the  boundary -layer  edge  must  be  established. 
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At  the  wall  we  apply  the  usual  no  slip  condition  and 


write 


u(^,(u,0)  - w(r,<jj,0)  = 0 


(16  ) 


To  allow  for  mass  injection  at  the  surface  we  write 

v(f,w,0)  » Vq(I,w)  (17) 

Using  the  subscript  m to  denote  matching  values,  we  write  for 
large  n, 

u--»Ujjj(C,w,n)  1 


w— Wjjj(C,aj,n) 

P-*-Pn,(C.(i).n) 


as  n — 


as ) 


In  order  to  find  values  for  u_  and  w_  we  turn  to  the 

m m 

requirement  of  vanishing  vorticity  for  large  n.  First  we 
note  that 


1 L.  ^ ~ 9V  3(h  w) 

n = V X V = < h^h  i. — 

“ " " h\h  1 ^ ^ 

C w I L 


3V  3(h.u;  ^ 3(h  w)  3(h_u)  1 

.hi  ~ ^ + i ? ^ ) 

“ “ L 3C  3n  J L 35  J J 


Substituting  Eqs.  (9)  and  (10),  the  boundary-layer  scaling 
laws  gives 
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Multiplying  by  e and  retaining  terms  to  0(e)  gives 


ij.  3(h  w)  i 3(h-u) 
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Requiring  0 to  vanish  at  large  n then  gives 
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Equations  (19 ) can  be  integrated  to  give 
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In  order  to  establish  Fj  and  Fg  the  matching  condition  that  the 
viscous  velocity  profile  as  n » match  the  inviscid  velocity 
profile  as  n — 0 will  be  used.-  Thus 


Fj^(C . w) 
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where  the  capitol  letters  denote  inviscid  quantities.  Ex~ 
pending  the  first  of  Eqs.  (21)  in  a MacLaurin  series  gives 
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The  outer  flow  must  be  irrotational,  however,  so  that  from  the 
definition  of  vortlcity 
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We  now  note  that,  by  a KacLaiirln  series. 
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Thus,  using  the  binomial  series 
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(1  - k^n  + k^n^)  - U(C,u>,0)(l  - k^n  + k^  n"")  + . . . 

*'c.0 

where  the  boundary  condition  V(\,ui,0)  * 0 has  been  used. 

Finally 

F^d.;;;)  - h^  q ua.w.o).  (22) 


In  a similar  manner 
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Equations (20)  then  become 
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since  h^  = h^  and  h = h . 

t,  t,  0)  (D 

3.  Elimination  of  Pressure  from  the  Boundary -Layer  Equations 
Examination  of  Eqs.  (11)  through  (15)  reveals  that 
the  principle  difference  between  these  equations  and  the  usual 
first-order  bound ary -layer  equations  is  the  fact  that  the  pres- 
sure is  not  constant  across  the  layer.  As  pointed  out  in  Ref. 

6 for  the  two-dimensional  case,  these  equations  can  be  made 
similar  to  the  first-order  boundary -layer  equations  by  elimina- 
ting the  pressure  as  a variable,  hence  eliminating  the  need  for 
direct  solution  of  the  normal  momentum  equation.  Once  this  is 
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accoiiq>li8hed,  solution  procedures  similar  to  those  for  the 
first -order  eqxiations  (see  Refs.  2 and  6 ) may  be  applied 
to  the  system. 


In  order  to  eliminate  the  pressxure  the  normal 
equation,  Eq.  (14)  is  first  integrated  to  give 
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where  use  has  been  made  of  the  matching  conditions  . 
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momentum 


(25) 

Now  note 

(26) 

(27) 


Eq.  (25)  is  differentiated  with  respect  to  w and  evaluated 
at  n-»-®  to  give 

8r  a[i^(0.;D  ^ i^o.n)]  ^ a RA\  , 

3ti)  3to  3u)  V Zhf  / 3(1)  V 2h  * 


3P 

+ — 
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In  order  to  evaluate  3p/9a)  as  n-*-®  the  n-^®  limit  of  Eq . (13)  is 
taken,  which  gives,  after  use  of  the  matching  and  inviscid  ir- 
rotationality  condition 

3a)  2 9(j,  \ / \ h^  / 

Use  of  this  result  in  Eq.  (28)  then  gives 
3F  3(l7+I-) 

3a)  3a) 


or,  upon  integration 
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Thus,  Eq.  (25)  becomes 
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In  a similar  manner,  this  equation  is  differentiated  v/ith 
respect  to  C and  evaluated  as  n — ® yielding 


II  = l£ 
H af 


(30) 


Evaluating  Eq . (12)  as  n — “ and  using  the  matching  and 
irrotationality  conditions  gives 
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Integration  of  this  result  with  respect  to  C and  substi- 
tution into  Eq.  (29)  gives 
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where  C is  a constant  of  integration.  In  order  to  evaluate 
C Eq . (31)  is  examined  as  n-<-<», 


p„(e.“> 


r2  ,,2  r ^2 

g , 0 e _ gi, 0 e 


(32) 


2h2 
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From  inviscid  conditions  near  the  body, 

W? 


P = ^+C 

e 2 2 


which  is  recognized  as  Bernoulli's  equation  and  hence  C is 
identified  as  the  stagnation  pressure.  Solving  Eq.  (32)  for 
C and  substituting  into  Eq . (31)  yields 


p(^,g),n)  = Pg(£;,g))  + — |1  - 


^ I + I-^(n,“>)  + I-(n,“)  (33) 

0)  ' 

The  pressure  is  finally  eliminated  from  Eqs.  (12)  and 
(13)  by  substitution  of  Eq . (33)  which  gives 
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4.  Transformation  of  the  Boundary-Layer  Equations 

We  now  define  the  following  new  independent  variables 

^ = I j 

0)  = 0) 


n(^,a),n)  = 

and  new  dependent  variables 
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where  W can  be  either  or  U . From  the  chain  rule 
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The  time  averaged  continuity  equation,  Eq . (11),  is  now 
to  give 


1 9 
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or,  after  using  Eqs.  (37)  through  (39) 
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This  equation  is  now  differentiated  with  respect  to  n 
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which  is  the  transformed  continuity  equation. 

The  convective  operator  is  now  transformed. 


u 9 w 9 9 h U F 3 ^ qU^F  9 

— — — z'*'  ''“z~  ■ ■ ' ^ — - — + 2 — 

h_  9^  h 9a)  9n  hr  9C  hr 

C 0)  i c. 


+ 


_WG  d h 'tIG  d I 2U  9 

“I — * -^^4-  "z  — *'  — - 

1 9a)  h 3n  \ C 9h 

0)  0) 


"I 


integrated 
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to  give 
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Using  Eq . (40)  this  result  becomes 

u9  w3  8 h^«UF9  h „WG  9 2U  9 

— — + — — + v — = + + _±  V — (43) 

h.  9C  h 9a)  9n  h;  9S  h 9a)  C 

The  transformation  relations,  Eqs . (37),  (38),  (39)  and  (43) 
are  now  applied  to  the  time  averaged  i momentum  equation,  Eq. 
(34)  yielding 
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The  transformation  relations,  Eqs . ( 37),  ( 38),  ( 39)  and 
( A2)  when  applied  to  the  time  averaged  to-momentimi  equation, 
Eq.  (35),  give 
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Finally,  the  transformed  time  averaged  energy  equation  is 
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5.  Transformed  Governing  Equations  in  the  Plane  of  Symmetry 


The  governing  equations  at  the  windward  symmetry  plane 
are  found  by  examining  Eqs . ( 37)  - (56)  in  the  limit  as  w— »0, 
In  the  plane  of  symmetry  we  take  W = and  note  that 


G = lim  — = lim 
to — — 0 0) — »"0 


(57) 


with  the  last  limit  being  bounded.  In  the  plane  of  symmetry 
we  note  the  following, 
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Substitution  of  Eqs.  (58)  to  (62)  into  Eq.  (42)  gives 
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which  is  the  windward  streamline  form  of  the  continuity 
equation.  We  further  note  that  at  the  symmetry  plane 

■ 0,  so  that  the  C-momentum  equation,  Eq.  (44),  becomes 
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Before  examining  the  o) -moment vim  equation  we  must  first 
look  at  some  limits  at  the  windward  streamline. 
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After  noting  that  both  of  these  quanitities  have  bounded  values ^ 
we  find  for  the  w-momentum  equat  ion,  Eq.  (52), 
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Now,  at  the  windward  streamline, 


9n  h 
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We  also  note  that 
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These  latter  results  allow  us  to  evaluate  the  last  term  in 
the  (j-momentum  equation. 

Finally,  we  write  for  the  symmetry  plane  energy  equation 
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V^e  also  note  that  at  locations  where  * 0 the  limit 

du) 

(57)  is  unbounded.  From  Eq.  (63)  we  also  see  the  term  requiring 

G is  not  required  and  hence  we  do  not  need  to  solve  the  w-momentum 

w 3W 

equation.  Also  note  that  ■ w and  hence  the  product 
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G — - is  bounded  at  the  symmetry  plane  when  - 0. 
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6.  Transformed  Equations  at  the  Stagnation  Point 


At  the  stagnation  point,  the  properties  at  the  edge 
of  the  boundary  layer  and  the  metrics  at  the  surface  are 
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Let  us  first  substitute  the  above  relations  into  the  inviscid 
irrotational ity  condition.  This  gives 
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We  now  substitute  these  results  into  the  continuity 
equation,  Eq . (42) 
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Thus , 
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Using  these  results.  Eq.  (44).  the  C -momentum  equation. 
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Before  the  w-momentum  equation  is  evaluated  note 
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With  these  results  the  w -momentum  equation,  Eq . (52) 


becomes 
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Finally,  the  energy  equation  becomes 
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7.  Coordinate  System  and  Metric  Coefficients 
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The  orthoponal  surface  curvilinear  coordinate  system 
used  to  write  the  boundary- layer  equations  is  identical  to 
the  one  developed  by  Blottner  and  Ellis  (Ref.  2).  As  shown 
in  Fie.  1 the  three  coordinate  directions  are  w,  n.  For 
each  set  of  oi,  n there  corresponds  a set  x,  6,  r in  the 
polar  coordinate  system  whose  origin  is  at  the  stagnation 
point.  The  inviscid  data  tape  contains  data  in  the  oj.n 
coordinate  system  and  at  each  point  the  corresponding  values 
of  X,  0,  r are  written.  The  transformation  between  the  two 
coordinate  systems  is  generated  numerically  in  program  BLOT 
which  is  contained  in  the  TAPCEN  program  (see  Ref.  7).  The 
procedure  is  identical  to  that  of  Ref.  2.  For  each  point  in 
the  0),  n system,  program  BLOT  also  generates  values  of 
h.  « and  h « and  these  quantities,  along  with  their  C"deriva- 

Cf  f \J  U) , u 

tives  are  also  witten  on  the  inviscid  data  tape. 

The  boundary-layer  program  itself  contains  routines 
to  generate  values  of  h^  and  h^  as  functions  of  n for  points 
in  the  boundary  layer  away  from  the  wall.  These  routines  are 
used  only  in  calculations  where  surface  curvature  effects  are 
included.  For  cases  where  no  surface  curvature  effects  are 
included,  the  following  relations  are  used: 
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In  order  to  calculate  the  values  of  the  metric 
coefficients  away  from  the  wall  the  definitions  of  the  metrics 
are  used,  i.e. 
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dC 
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where  s and  t are  defined  in  Fig.  2 The  following  definitions 
are  now  made , 
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The  above  expressions  are  evaluated  numerically 
using  the  values  of  x,  r,  and  <[)  for  points  in  the  profiles. 
If  g is  taken  to  represent  x,  r or  (p  then  (see  Fig.  3). 
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The  positions  of  the  coordinates  themselves  are 
generated  with  the  numerical  scheme  of  Blottner  and  Ellis 
(Ref.  2 ).  This  procedure  is  first-order  accurate  and  re- 
quires relatively  small  step  sizes  in  ^ and  w.  For  this  reason 
the  coordinate  generation  routine  was  coded  independently  of 
the  boundary- layer  program  so  that  the  coordinate  system  could 
be  generated  on  the  fine  mesh  and  the  boundary -layer  equations 
evaluated  over  a coarser  mesh. 

8.  Eddy  Viscosity  Models 

Prandtl's  mixing  length  hypothesis  states  that  the 
eddy  viscosity  is  the  product  of  some  characteristic  length 
and  the  normal  velocity  gradient.  The  characteristic  length 
is  related  to  the  size  of  the  eddies  of  momentum  flux  normal 
to  the  body  and  is  called  the  mixing  length.  For  two-dimensional 
flow  this  concept  leads  to: 

e •=  p^*^|3u/3n|  (88) 
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Prandtl's  studies  assumed  that  the  eddy  viscosity  should 
depend  only  on  local  eddy  scale  and  on  the  properties  of  tur- 
bulence. Adams  (Ref.  8 ) extended  this  concept  to  the  three- 
dimensional  case  by  assuming  that  the  eddy  viscosity  is  also 
independent  of  coordinate  direction  by  writing  the  component 
of  turbulent  stress  terms  as: 

T = -pu’v'  = 9E/3n  9u/9n  (89) 


Tj.  = -pv'w'  = olj-  8E/Sn  3u/3n  (90) 

(t) 

where  E is  some  scalar  function.  Therefore, 

e = €-  = e = pS.*^  3E/3n  (91) 
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The  total  shear  in  each  direction  is  written  as: 

= y 9u/3n  - pu'v'  = y 3u/3n  + 3u/3n  (92) 

T = y 3w/3n  - pw'v'  = y 3w/3n  + e,,  3w/3n  (93) 

CO 


therefore  the  total  resultant  shear  is  written  as: 
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Using  equations  (94)  and  (91)  the  total  resultant 
shear  becomes: 

T » j^y  + pi*^  3E/3nj  [(3u/3n)^  +(3w/3n)^  j ^ (95) 
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By  analogy  with  the  two-dimensional  case  where  the  eddy 
viscosity  expression  incorporates  the  velocity  gradient  of 
the  shear  component,  the  scalar  E becomes: 


3E/3n  = 


^3u/9n)^  +(3w/3n)^J 


1/2 


(96) 


0) 
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which  reduces  to  the  two-dimensional  form  when  w “ 0.  This  is 
referred  to  as  the  invariant  turbulence  model  by  Hunt,  Bushnell, 
and  Beckwith  (Ref.  9 ),  and  was  used  with  success  by  Adams 
(Ref.  8 ). 

The  model  used  in  this  investigation  is  the  common  two- 
layer  inner-outer  model  which  uses  the  Prandtl  mixing  length| 
theory  and  the  Van  Driest  or  Peichardt  damping  near  the  wall. 
Following  Patankar  and  Spalding  (Ref.  10  ) and  Adams  (Ref.  8 ) 

the  mixing  length  distribution  is  as  follows: 

n {0  < n ^ Arij^/k^} 


{Anj^/k*  < n} 


(98  ) 


where 


k*  - 0.435 
\ - 0.09 

- n when  [^(u^  + w^)y(ug^  + W^^)! 


0.99 


The  inner  law  is  damped  near  the  wall  so  as  to  yield  the  exact 
laminar  shear  stress  term  at  the  wall.  To  accomplish  this,  two 
different  damping  factors  have  been  used  in  this  investigation. 
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Van  Driest' s damping  term  with  local  shear  stress,  and 
Reichardt's  (Ref. 11  ) damping  term. 


is : 


Van  Driest's  damping  term  for  two-dimensional  flow 


(99) 


where  x is  the  local  shear  stress  and  A is  26.0.  Therefore 
the  total  shear  near  the  wall  becomes: 


2 2 

T “ ij  3u/3n  + pk^  n 


1 - exp 


(3u/3n)  (100) 


for  two-dimensional  flow.  Again,  use  is  made  of  analogy  to 
derive  the  form  of  the  near  wall  shear  for  a three-dimensional 
flow.  By  analogy  of  equation  (100)  with  equations  (95)  and 
(96)  the  three-dimensional  form  of  the  total  shear  becomes 

2 


or 


2 2 

= p 9E/9n  + pk*  n 


- I,  2 2 

- pk^  n 


1 - exp 


^jWTp 

yA 


(9E/9n)  (101) 


(9E/9n) 


(102) 


Cebeci  (Ref.  12)  developed  a mass -transfer  correction  to 

Van  Driest's  inner  eddy  viscosity  law  by  modifying  the  damping 
★ 

constant  A . For  turbulent  flows  with  mass  transfer  Cebeci 
determined  the  damping  constant  to  be 

A*  - 26  exp  (-5.9  Vq  "^) 
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where 


Reichardt ' s expression  for  the  inner  eddy  viscosity 
law  was  obtained  by  curve  fitting  experimental  pipe  flow 
data.  The  expression  is: 

. uk*  . 11.0  tanh  ('^^  )]  (103) 

As  can  be  seen  this  expression  does  not  involve  the  velocity 
gradient  terms.  For  this  reason  it  is  preferred  for  use  in 
numerical  solutions,  since  it  usually  requires  fewer  itera- 
tions to  converge. 

Following  equations  (97)  and  (98)  the  outer  eddy 
viscosity  law  is: 

Sq  = 9E/:n  (104) 


and  the  total  shear  stress  is: 

Tq  - M 3E/9n  + X^  (9E/9n)^ 


(105) 


The  outer  eddy  viscosity  law  is  used  in  conjunction  with  the 
Klebanoff  (Ref.  13)  intermittency  factor  which  assures  a 
smooth  approach  of  Cq  to  zero  as  y -*•  6.  The  modified  law  is: 

Cq  - X^  y 9E/9n  (106) 


where  y is  Klebanoff 's  interm 
Y - [l  + 5.5  (n/6)^ 


ttency  factor: 
-1 


(107) 
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Schetz  and  Favin  (Ref. 14  ) have  derived  a correction  to 
Reichardt's  inner  eddy  viscosity  law  for  cases  of  mass  trans- 
fer. This  correction  has  been  used  in  the  current  investiga- 
tion, giving  this  corrected  expression  for  the  inner  eddy 


viscosity: 


= ku  (1  + Vq^  (n*"  - n^"*"  tanh  (n'^/n^'*') ) (108) 


/ « +\ 


where 


n^  = n^Tp/\i 


n^"^  = 3.65/(Vq'*'  + 0.344) 


The  quantity  u is  found  by  integration  of  the  expression 


(1  + v^'*’  u'*') 


1 + k (1  + Vq^  u^)  (n"*"  - n^^  tanh  (n'^/n^'^)) 


(109) 


or  using  equation  (108) 


du"^  = + V u'*') 

dn^  (1  + e^) 


(110) 


Since  the  eddy  viscosity  is  implicit  in  the  integration  for 
u^,  the  calculation  of  is  an  iterative  procedure  for  mass 
transfer  cases. 
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9.  Transition  Models 


I 

f 


f 


Two  models  of  transition  from  laminar  to  turbulent  flow 
have  been  used  in  this  investigation.  One  model  is  a simply 
instantaneous  transition  to  turbulent  flow,  and  there  really 
is  no  transition  region  or  zone  at  all.  In  the  second  case  a 
smooth  transition  to  turbulent  flow  occurs  over  a prescribed 
distance.  This  distance  is  known  as  the  transition  zone 
and  is  defined  as  the  distance  between  the  onset  of  transition 
at  C “ and  the  beginning  of  fully  turbulent  flow  at  C * Cj. 
at  some  point  downstream. 

The  probability  of  turbulent  flow  at  any  point  is  ex- 
pressed by  a model  by  Dhawan  and  Narasimha  (Ref . 15  ) as: 

= 1 - exp  (-  4>  ((X-X^)/x)^)  (111) 

where  I^Cx)  is  the  transition  intermittency  factor, 
and 

$ - 0.412 

X - Xi  - Xj 

4 = 0.75  ^f  - 0.25 


and  where 


If(X^)  - 0 

Ij(X^)  - 0.97 


(112) 


By  substituting  equation  (112)  into  (111)  an  expression  for  x can 
be  found  based  on  the  transition  zone  length: 

7 “ (X.J.  - Xj.)/2.917  (113) 
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Now,  substituting (113)  back  into  (111)  the  final  expression 
for  the  transition  intermittency  factor  as  used  in  this  in- 
vestigation is  obtained: 

lf(0  = 1 - exp  [ 0.412  (2.917)^  ] (114) 

The  transition  intermittency  factor  is  employed  as  a simple 
multiplier  of  the  eddy  viscosity  in  the  governing  equations 
and  therefore  acts  as  a damping  coefficient  for  the  full  tur- 
tulent  eddy  viscosity.  It  is  an  expression  relating  the 
fraction  of  time  any  particular  point  spends  in  turbulent  flow, 
and  therefore  the  probability  of  turbulent  flow  existing  at 
that  point. 
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10.  Finite-Difference  Method 

The  finite  difference  method  used  in  this  investiga- 
tion is  idential  to  the  one  used  by  Frieders  and  Lewis  (Ref. 17  ). 
Basically  the  procedxrre  is  based  on  the  method  of  Dwyer 
(Ref . 18  ) with  modifications  by  Krause  (Ref. 19  ),  The  pro- 
cedure allows  for  variable  spacing  of  the  normal  coordinate. 

As  implemented  in  the  current  investigation,  the 
finite-difference  procedure  is  a forward  marching  one,  thus 
taking  advantage  of  the  parabolic  nature  of  the  governing 
equations.  The  method  used  marches  away  from  the  stagnation 
point  by  first  stepping  down  the  windward  symmetry  plane  one 
step  and  then  marching  around  the  body  from  the  windward  to 
leeward  symmetry  plane.  This  process  is  repeated  until  the 
calculation  is  completed. 

Basically  the  finite-difference  procedure  is  implicit 
in  the  normal  direction  and  explicit  in  the  and  co-directions. 

In  order  to  retain  stability  in  regions  of  reversed  cross -flow, 
the  Krause  difference  molecule  is  used  for  w differencing  as 
depicted  in  Fig.  3 . Taking  w to  be  a general  variable  (i.e. 
either  F,  G or  6)  then 


9w 


w 


2^ 


w 


3 .n 


A5 


^ . <“2.n  - "l.n>  <"4.n  ' '*3.n> 

3co  2Lbi 
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.V 


where  the  subscript  notation  is  that  of  Fig.  3 . Central 
differences  are  used  in  the  normal  direction  with  mesh  points 
spaced  according  to  the  formula 


k = 


An 


n 


An. 


n-l 


where  k is  a constant  which  can  be  set  at  the  discretion  of 
the  user.  Substitution  of  the  finite-difference  expressions 
into  the  governing  equations  results  in  a set  of  non-linear 
difference  equations  of  the  form 


■V2 


- 4“  B \jjn 
,n-l  n'^2,n 


^ 2 - 
-C  Wo  ,1  +E  Wo  = D 

n 2 ,n+l  n 2 ,n  n 


(115) 


This  relation  is  linearized  using  the  Newton-Raphson  interation 
formula 


2 0 0 / 0 v2 

w.,  ^ = 2W9  „ w,  ^ - (w,  ) 
4.n  4 ,n  4 ,n  4 ,n 


(116) 


where  w^  ^ is  the  value  of  the  dependent  variable  from  the 
previous  iteration.  For  the  initial  iteration,  w^  ^ is  ap- 
proximated v/ith  V7,  . Use  of  Eq.(116)  in  Eq.(115)  results  in 

1. , n 

a set  of  simultaneous  linear  algebraic  equations  of  the  form 


■A  w , +B  w„  -M2  w.,  . 1 = 

n 2, n-l  n 2,n  n 2,n-fl  n 


which  are  solved  using  the  Thomas  algorithm  (Ref.  20). 
11.  Normal  Pressure  Gradient  Approximation 

The  normal  pressure  gradie  it  which  appears  in  the 
normal  momentum  equation,  Eq.  (14)  basically  arises  from 
centrifugal  force  effects.  Due  to  the  low  speeds  being 
considered  in  the  current  work  this  effect  was  neglected. 
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This  was  done  by  setting  and  to  zero  in  Eqs.  (44) 
and  (52) . In  the  actual  coding  however  these  terms  were 
left  in  the  governing  equations  and  this  effect  could  easily 
be  included  in  the  future  by  having  subroutine  PRESSI  evaluate 
Eqs.  (45)  and  (46)  in  program  ICBL3D. 


Section  III 


RESULTS  AND  DISCUSSION 

In  order  to  test  the  validity  and  accuracy  of  the 
computer  program  developed  from  the  foregoing  analysis,  a 
variety  of  test  cases  was  selected.  Since  the  principle  de- 
parture of  the  present  analysis  from  previous  efforts  is  in 
the  inclusion  of  surface  curvature  effects,  it  was  this  area 
that  received  considerable  attention  during  the  testing  process. 
Test  cases  were  run  at  a Reynolds  number  of  one  hundred  in 
order  to  achieve  a thick  boundary  layer  on  the  body,  thus 
amplifying  the  curvature  effects.  Calculations  were  made  over 
spheres  and  spheroids  at  various  angles  of  attack  and  turbulent 
effects  were  included  in  one  comparison. 

The  first  test  case  run  was  that  of  a sphere  at  angle 

of  attack.  In  this  calculation  a coordinate  system  is  selected 

that  is  not  aligned  with  the  sphere's  wind  axis  thus  producing 

the  requirement  for  a three-dimensional  calculation  of  the  body 

boundary  layer.  The  utility  of  this  calculation  is  that  it  can 

be  easily  compared  with  existing  axisymmetric  boundary-layer 

calculations.  Calculations  were  made  for  a unit  sphere  at 

a = 2°,  Re  = 100.  The  results  were  compared  with  those  of 
00 

Davis  et  al.  (Ref.  6).  The  results  of  Ref.  6 were  obtained 
with  an  axisymmetric  boundary -layer  code  with  longitudinal  and 
I transverse  curvature  effects  included  (SFC) . The  results  of  these 

' comparisons  are  shown  in  Figs.  4 through  6.  Plotted  in  Fig.  4 

is  the  development  of  the  skin  friction  and  displacement  thickness 
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along  the  windward  S3nnmetry  plane  of  the  sphere.  The  results 
of  the  present  calculation  without  surface  curvature  (NoSFC) 
agree  identically  with  predictions  made  by  the  VPI61SU  axisymmetric 
boundary- layer  code  (Ref.  16).  The  skin- friction  comparison 
between  the  present  method  and  that  of  Davis  et  al.  (Ref.  6) 
for  the  SFC  case  is  also  excellent.  Figure  5 shows  the  variation 
of  skin  friction  around  the  sphere  along  ^“constant  lines.  Com- 
pared in  Fig.  6 are  the  velocity  profiles  at  a point  on  the 
windward  symmetry  plane  for  the  with  and  without  SFC  cases. 

Calculations  were  then  made  over  a 4:1  spheroid  at 
a = 2°  and  Re^  = 100.  The  semi -major  axis  of  the  body  was 
aligned  with  the  angle  of  attack  line  and  was  four  feet  long. 

The  results  of  these  calculations  are  shown  in  Figs.  7 through 
10.  Figure  7 illustrates  the  development  of  the  boundary  layer 
along  the  windward  symmetry  plane.  The  skin  friction  for  the  No 
SFC  case  is  compared  with  an  unpublished  calculation  made  with 
the  Blottner  and  Ellis  (Ref.  2)  code.  As  can  be  seen  the  agree- 
ment is  excellent.  Figure  8 illustrates  the  development  of  the 
skin  friction  around  the  body  and  compares  the  results  with 
that  of  the  Blottner  and  Ellis  code  for  the  NoSFC  case.  In 
Figs.  9 and  10  velocity  profiles  are  compared  between  the  two 
codes . 

An  interesting  result  of  these  calculations  is  that 
inclusion  of  surface  curvature  effects  initially  produces  a 
somewhat  thicker  body  boundary  layer  compared  with  the  NoSFC 
case.  This  is  in  contrast  to  the  inclusion  of  transverse  curva- 
ture effects  only  which  initially  produces  a thinner  boundary 
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layer  compared  with  the  NoSFC  case  (see  Ref.  6).  Since  the 
pressure  gradients  experienced  by  the  boundary  layer  over 
flatter  portions  of  the  spheroidal  body  where  longitudinal 
curvature  effects  become  less  important  are  not  such  as  to 
thin  the  boundary  layer,  this  thicker  boundary-layer  persists 
along  the  body.  The  results  indicate  that  all  curvature  effects 
must  be  accounted  for  in  order  to  correctly  predict  the  boundary 
layer  growth  over  blunt  nosed  bodies. 

Figures  11  and  12  detail  the  boundary -layer  growth 
for  the  same  4:1  spheroidal  body  at  a = 10°.  The  effects  of 
this  larger  angle  of  attack  are  clearly  seen  by  comparing  the 
circumferential  skin  friction  plot  of  Fig.  12. 

In  Figure  13  a turbulent  calculation  is  presented  for 
a 4:1  spheroid  at  a = 0°  and  Re^  = 10^.  Comparison  is  made 
between  the  present  method  and  the  results  of  Chang  and  Patel 
(Ref.  3).  As  can  be  seen  the  agreement  is  excellent. 

These  calculations  were  made  to  verify  the  code  and 
as  can  be  seen  by  the  calculated  results,  agreement  between 
the  present  method  and  previously  published  works  is  excellent. 
It  should  be  borne  in  mind , however , that  the  present  code  is 
much  more  versatile  than  either  the  codes  of  Refs.  2 or  3 in 
that  it  includes  surface  curvature  effects,  it  can  treat 
arbitrary  body  shapes  and  it  includes  turbulence  effects. 
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SECTION  IV 


r 


CONCLUDING  REMARKS 

The  total  system  of  computer  programs  generated  to 
implement  the  foregoing  analysis  represents  a package  with 
very  general  capabilities.  These  capabilities  are  as  follows: 

1.  Body  must  have  a blunt  nose  and  a plane  of  symmetry 
but  is  otherwise  arbitrary  so  far  as  the  boundary- 
layer  code  is  concerned  at  arbitrary  angle  of  at- 
tack . 

2.  Current  inviscid  capabilities  are  restricted  to 
axisymmetric  shapes  at  zero  and  non-zero  angles 
of  attack  and  arbitrary  cross-sections  at  zero 
angle  of  attack. 

3.  All  surface  curvature  effects  included. 

4.  Laminar,  transitional  and/or  turbulent  flows  can 
be  calculated. 

5.  Effects  of  heat  and  mass  transfer  included. 

To  the  authors'  knowledge  these  capabilities  represent 
the  most  complete  package  available  today  for  predicting  three- 
dimensional  boundary  layers. 

The  system  of  programs  was  separated  into  two  independent 
systems  of  programs  in  order  to  maximize  operational  versatility 
and  to  facilitate  future  development.  By  removing  the  calcu- 
lation of  the  inviscid  flow  and  coordinate  system  from  the 
boundary-layer  code,  this  boundary- layer  code  can  be  viewed  as 
a solution  procedure  for  the  set  of  governing  equations 
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developed  in  Section  II.  The  separate  inviscid  flow  computer 
program  supplies  data  to  this  solver  which  defines  both  the 
body  geometry  and  inviscid  flow.  Thus,  if  alternate  means  for 
generating  either  the  inviscid  flow  or  coordinate  system  are 
desired,  other  methods  can  be  substituted  into  the  appropriate 
blocks  in  the  inviscid  package  without  affecting  the  boundary- 
layer  code.  Further,  the  separate  boundary- layer  calculation  is 
allowed  to  proceed  on  its'  own  step  size  along  the  surface  with- 
out considering  mesh  requirements  of  the  inviscid  flow  or  co- 
ordinate generation  codes. 

In  summary,  the  analysis  and  programs  resulting  from 
this  investigation  represent  as  versatile,  flexible,  general 
and  efficient  a method  for  predicting  incompressible  three- 
dimensional  boundary-layers  available  today. 
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Figure  4.  BOUNDARY- LAYER  DEVELOPMENT  ALONG  WINDWARD 
STREAMLINE  OF  A SPHERE  AT  a = 2° 
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Figure  5.  SKIN  FRICTION  DISTRIBUTION  AROUND 
A SPHERE  AT  a = 2° 
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Figure  9.  COMPARISON  OF  BOUNDARY  LAYER  VELOCITY 
PROFILE  ON  4:1  SPHEROID  at  a - 2° 


Figure  10.  COMPARISON  OF  BOUNDARY -LAYER  VELOCITY 
PROFILES  4:1  SPHEROID  at  a - 2° 
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Figure  12.  SKIN  FRICTION  AROUND  A 4:1  SPHEROID 
at  It  * 10° 
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TURBULENT  BOUNDARY -LAYER  ON  A;1  SPHEROID  AT  ZERO  ANGLE  OF  ATTACK 


